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Planetesimal formation via fragmentation in 
self-gravitating protoplanetary discs 
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ABSTRACT 

An unsolved issue in the standard core accretion model for gaseous planet formation is 
how kilometre-sized planetesimals form from, initially, micron-sized dust grains. Solid 
growth beyond metre sizes can be difficult both because the sticking efficiency becomes 
very small, and because these particles should rapidly migrate into the central star. 
We consider here how metre sized particles evolve in self-gravitating accretion discs 
using simulations in which the gravitational influence of the solid particles is also 
included. Metre-sized particles become strongly concentrated in the spiral structures 
present in the disc and, if the solid to gas density ratio is sufficiently high, can fragment 
due to their own self-gravity to form planetesimals directly. This result suggests that 
planetesimal formation may occur very early in the star formation process while discs 
are still massive enough to be self-gravitating. The dependence of this process on the 
surface density of the solids is also consistent with the observation that extrasolar 
planets are preferentially found around high metallicity stars. 
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1 INTRODUCTION 

There are currently two models for the formation of gaseous, 
Jupiter- like plane ts. The most widel y accepted is the core 
accretion model (iPollack et alJll99d) in which, initially, a 
core of rock/i ce grows via th e collisional accumulation of 
planetesimals iSafronovl Il972l) and then, once sufficiently 
massive, it accrete s a gaseous envelope llLissauerl Il993t 
iPollack et aflll996h . The alternative scenario does not re- 
quire the initial growth of a core and assumes that pro- 
toplanetary discs may become gravitationally unstable and 
that gas giant planet s may subsequently form via direc t 
gravitational collapse (lBosslll99gl . 120001 iMaver et alJl2004) . 

One aspect of the core accretion model that has yet to 
be satisfactorily resolved is how the kilometre-sized planetes- 
imals, that ultimately coagulate to form the plane tary core, 
grow from, initially, micron-sized du st grains (e.g. ISafronovl 
Il972t IWeidenschilling fc Cuzzilll993h . Core accretion mod- 
els generally start with the assumption that these planetes- 
imal have already formed. However, solid growth beyond 
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metre sizes can be very difficult, due to two effects. On the 
one hand the sticking efficien cy of solids becomes relatively 
small in this size range (e.g- ISupulver et aljfl997h . On the 
other hand, it is well known that solid particles are influ- 
enced by gas drag and, for standard disc geometries, will 
generally lose angular momentum and migrate in towards 
the central star at a r ate that depends on the particle size 
iWeidenschillinell977f) . For a circumstellar disc with proper- 
ties appropriate for planet formation, the maximum inward 
radial velocity will generally occur for particles with sizes 
between 1 cm and 1 m and may easily exceed 10 3 cm s _1 
JWeidenschillindll977lh With such large inward radial ve- 
locities, it is possible that these objects may drift into the 
central star before becoming large enough to decouple from 
the disc gas, preventing the growth of the planetesimals re- 
quired for the formation of the planetary cores. 

The instability model, since it does not require the for- 
mation of a core, is unaffected by the above process. This 
scenario does, however, require that the disc be relatively 
massive, an d that the disc be able to coo l extremely effi- 
ciently (e.g. lGammie!l200il : lRice et alJl2003ft . Although very 
few Class II protostars appear to hav e discs with the re- 
quired mass (IBeckwith &c Sargenfll99ih . recent observations 
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su ggest that massive disc s may be prese nt during the Class 
^Rodriguez et alJl2005l) and Class I (lEisner et alJ l2005l) 
phases. Although this suggests that massive discs may be 
present at some stage during the star formation process, it 
is still not clear that such discs can cool sufficiently rapidly 
for the direct formation of gas giant planets. 

Even if gas giant planets cannot form via direct grav- 
itational collapse in discs around Class and Class I pro- 
tostars, it is still likely that these discs will experience self- 
gravitating phases i n which spiral structures will develop 
l)Lin fc Pringlelli98^ . The gas pressure on the inner edges 
of these spiral waves increases with radius, leading to super- 
Keplerian gas velocity. The drag force between the solid 
particles and the gas then results in the dust grains and 
small particles drifting towards th e peaks of the spiral struc- 
tures ( Haghig hipour fc Bossl2003t) . rather than simply drift- 
ing towards the central star. lRice et alJ i2004h have already 
shown, using numerical simulations, that for certain parti- 
cles sizes, the local density could be enhanced by a factor 
of 10 or more. A similar eff ect would occur in the presence 
of rings (D urisen et al.|2005l) . vortices jGodon fc Liviol2000t 
iKlahr fc Bodenheimerl2003f) and it has also been shown that 
density variations in the disc gas resulting from magneto- 
rotational turbulence can a lso produce significan t enhance- 
ments in particle densities llJohansen et al 

Particle density enhancements may aid planetesimal 
formation in two ways. Firstly, the enhanced collision rate 
l|Rice et alJ 1200 4) could aid collisional grain growth, and 
secondly, the densities achieved may be sufficient for plan- 
etesimal formation through direct gravitational collapse of 
the solid component of the disc (Goldrcich & Ward 1973; 
lYoudin fc Shull2002l) . It is this second possibili ty that we 
invest igate in this paper. We extend the work of iRice et alJ 
( 2004) to include the gravitational influence of the solid par- 
ticles. We consider three-dimensional, global simulations of 
self-gravitating accretion discs in which the gas is main- 
tained in a state of marginal gravitational instability and 
is coupled to the solid particles via a drag force. We include 
the gravitational influence of both the gas and the solid par- 
ticles self-consistently and we consider different particle sizes 
and different initial gas to solid particle surface density ra- 
tios. We find that if the initial particle surface density is 
sufficiently high, a significant fraction of the solid compo- 
nent of the disc, which is strongly concentrated in the self- 
gravitating spiral waves, can become unstable and subse- 
quently fragment into bound objects. We therefore conclude 
that what used to be considered the most difficult step in 
planetesimal formation (i.e. the growth beyond metre sizes) 
could actually be very rapid, thus leading directly to kilo- 
metre sized objects. 



2 NUMERICAL SIMULATIONS 

T he simulations per formed here are very similar to those 
of IRice et alJ |2004). We consider a system comprising a 
point mass, representing the central star, surrounded by two 
interpenetrating discs, a gas disc, and a 'planetesimal' disc. 



2.1 The gas disc 

The three-dimensional gaseous disc is modelled using 
smoothed particle hydrodynamics (SPH), a La grangian hy- 
drodynamics code (Bcnz 1990; iMonagh an 1992), and is rep- 
resented by 250000 pseudo-particles. Each particle has a 
mass, an internal energy and a smoothing length that varies 
with time to ensure that the number of neighbours (SPH 
particles within two smoothing lenghts) remains ~ 50. These 
neighbouring particles are used to determine the gas density 
which, together with the internal energy, is used to deter- 
mine the gas pressure. A tree is used to determine gravita- 
tional forces, and to determine the gas particle neighbours. 

The simulation starts by considering only the gaseous 
disc which has a mass Mdi sc =0.25 Mq and extends from 
0.25 to 25 au around a 1 Mq central star. The initial surface 
density profile R 1 , the initial temperature profile 

is T oc R~ ' 5 , and the temperature is normalised s uch that 
the d isc is initially gravitationally stable at all radii llToomrel 
1964). The disc gas is allowed to heat up through both PdV 
work and viscous dissipat ion with the visco sity given by the 
standard SPH viscosity jMonagha nl ll992ll . In the absence 
of cooling the gas has an adiabatic equation of state with 
7 = 5/3. 

The disc gas is cooled with a radially dependent cool- 
ing time of t C ooi = /3f2 -1 where Q, is the orbital angu- 
lar frequency, and /3 is a constant that determines the 
cooling rate relative to the angular frequency. In thermal 
equili brium, the cooling ti me and the viscous a parameter 
JShakura fc Sunvaevlll973T) . that characterises angular mo- 
mentum transport, are related through 
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where 7 is the ra tio of the specific heats. It has been shown 
jRice et alJl2005Tl that if the value of a required for thermal 
equilibrium is too large (a >~ 0.06) the disc will fragment 
rather than settle into a quasi-steady, long-lived state. We 
therefore use /3 = 7.5 which ensures that the related a value 
is small enough for the disc, in the simulations presented 
here, to settle into a long-lived, quasi-steady, self-gravitating 
state <Gammidl200lt IRice et alJl2003t lLodato fc Ric3l2004L 
l2005t IRice et alJl2005l) . The 'planetesimal' disc is added to 
the simulation once this long-lived state has been achieved. 



2.2 The 'planetesimal' disc 

The 'planetesimal' disc is also modelled using SPH and is 
represented by 125000 pseudo-particles. It is, however, as- 
sumed to be pressureless and hence the solid particles have 
no internal energy and experience only gravitational forces, 
and a drag force due to the differen ce between their veloc - 
ity and the velocity of th e disc gas llW eidenschilling 1977). 
As discussed in detail in IRice et alJ i2004f) . the drag force 
depends on the particle size, the local gas density, the local 
gas velocity, and the local gas sound speed. To determine 
the drag force each 'planetesimal' particle has a 'smoothing 
length' that is varied to ensure that the number of gas par- 
ticle neighbours remains ~ 50. These neighbouring particles 
are then used to calculate the gas density, velocity, and in- 
ternal energy at the location of ever y 'planetesimal' p article 
using the standard SPH formalism (lMonaghanlll992l) . The 
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exact value of the drag force is then determined by speci- 
fying the size of the solid particle. In each simulation, the 
'planetesimal' disc is assumed to contain particles of a single 
size. 

The primary differ ence between the work presented here 
and that presented bv lRice et al.l (120041) is that we include 
here the gravitational influence of the solid particles, rather 
than treating them as test particles. As with the gas parti- 
cles, the gravitational force due to the 'planetesimal' parti- 
cles is determined by including them i n the tree. T he details 
of the tree method can be found in iBenj il99(j) . but es- 
sentially distant particles are grouped into nodes and their 
gravitational force is computed from the position of the node 
and its quadrupole correction. Very nearby particles are in- 
cluded in the neighbour list of the particle being considered, 
and the force due to these neighbouring particles is calcu- 
lated directly. A direct gravitational force calculation that 
includes a 'planetesimal' particle is also softened. This is 
needed because 'planetesimal' particles do not feel any pres- 
sure forces and hence may undergo extremely close encoun- 
ters with other particles. In most of the simulations we used 
a softening of 10 -2 au, but repeated some simulations using 
1CT 1 au and 1CT 3 au. 

The 'planetesimal' disc is only introduced once the 
gas disc has settled into a long-lived, quasi-steady, self- 
gravitating state. In this self-regulated state, the gravita- 
tional stability parameter, Q ~ 1, and the disc temperature 
increases with radius from ~ 30 K in the inner disc to ~ 150 
K in the outer disc. The disc temperature is therefore well 
below the dust sublimation temperature of ~ 1600 K and 
we would not expect grain evaporation to be an important 
process in the self-gravitating regions of the disc. It should, 
however, be noted that the magnitude of the temperature 
and the temperature profile are largely determined by our 
choice of disc mass and surface density profile. However, it 
would require a disc almost two orders of magnitude more 
massive than the one in these simulations for grain evapo- 
ration to be an important effect. A slightly more massive 
disc could, however, have a temperature that may be high 
enough to influence the condensation of ices. 

The 'planetesimal' disc is initially located only in the 
midplane (z — 0) and extends from 2 to 20 au. The particles 
are distributed randomly in such a way as to give the initial 
surface density profile of E p i oc R~ , the same as that of the 
initial gas disc. An initial spiral structure is not introduced 
into the 'planetesimal' disc. Although the disc is initially 
located only in the midplane, it is given an isotropic velocity 
dispersion with a magnitude of 0.1c s where c s is the local 
gas sound speed. This ensures that the planetesimal disc is, 
at the start of the simulation, gravitationally stable. The 
'planetesimal' particles are free to move vertically allowing 
the disc to achieve a finite thickness, which turns out to be 
comparable to the gas disc thickness. 

Once the 'planetesimals' have been added, the simula- 
tion is evolved for about an additional outer rotation period 
(~ 125 yrs). We consider particle sizes of 150 cm and 1500 
cm, and gas to dust surface density ratios of 100 and 1000. 




Figure 1. Surface density structure of a 0.25 Mq disc around a 1 
Mq star that has reached a long-lived, self-gravitating state. The 
scale shows the logarithm of the surface density, E gas , with the 
scale covering 2.4 < log(E gas /gcm -2 ) < 4.4. 

3 SIMULATION RESULTS 

As discussed above, the simulation initially considers 
only the gaseous disc component. This is evolved until 
it settles into a m arginally stable, self-gravitating state 
jLodato fe Ricell2004D . In this quasi-steady state, the in- 
stability produces spiral-like structures in the disc. This 
is illustrated in Figure Q which shows the surface density 
structure of the gaseous disc, once it has reached this long- 
lived, self-gravitating state. Although it is not clear that 
a protostellar disc could be self-gravitating at radii of a 
few au, it is quite likely that it plays a vital role in trans- 
porting an gular momentum in the earliest phase s of star 
formation (|Lin fc Pringlelll990l) . feell fc Linl l|l994Fl also in- 
voke self-gravity at ~ 1 au in their models of FU Orionis 
outbursts, and the p resence of a dead zone iGammielll996l : 
lArmitaee et al]|200ll) could also allow a self-gravitating re- 
gion to extend to small radii. 

Once the gas disc has evolved into this self-gravitating 
state, we introduce the 'planetesimal' disc as discussed in 
For the gas disc we consider here, the largest radial drift 
rates will occur for particles with sizes of ~ 100 cm (see 
iRice et aU 12004ft . It is particles of this size that we expect 
to bec ome the most strongly concentrated in the spir al struc- 
tures jHaehighipour fc BosJl2003t IRice et al.ll2004ft . Parti- 
cles significantly larger, or significantly smaller, would not 
be as strongly influenced by the spiral structures. 

We consider initially 150 cm particles with a surface 
density 100 times smaller than that of the gas disc. Fig- 
ure |5| shows the surface density structure of this disc ~ 80 
years after these particles are introduced. The colour scale 
differs from that in Figure Q due to the different initial sur- 
face densities. What is clear from Figure|5|is that the spiral 
structures have fragmented into clumps. These clumps are 
bound and have densities significantly larger than the that 
of the disc gas. We also find fragmentation throughout the 
planetesimal disc. At the time shown in Figure^ the inner- 
most clump is at r = 2.1 and the outermost is at r = 16.2. It 
is likely that the outermost regions of the 'planetesimal' disc 
would also produce clumps if we were to run the simulation 
further. 

The clump formation starts extremely quickly. Figure[3] 
shows the fraction of solid particles that have a local density 
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Figure 2. Surface density structure of the 'planetesimal' disc 
consisting of 150 cm particles ~ 80 years after the 'planetes- 
imal' disc is introduced into the simulation. The scale shows 
the logarithm of the surface density, £ p i, with the scale cover- 
ing —0.6 < log(£ p i/gcm~ 2 ) < 2.4. In this simulation, the spiral 
structures in the 'planetesimal' disc have fragmented producing 
bound clumps. 



exceeding 6 x 10~ 8 g cm -3 . This threshold density exceeds 
the gas density everywhere in the disc and so ensures that 
these particles are all located in dense clumps. The first ev- 
idence of clump formation occurs ~ 20 years after the solid 
particles are introduced into the simulation, and almost 15% 
of the 'planetesimal' particles are located in dense clumps 
after ~ 80 years. Although this gives some idea of the effi- 
ciency of the process, we were unable to continue the simu- 
lation much further due to the exceedingly small timesteps 
at which particles in dense clumps have to be advanced. 
We therefore cannot assess if the clump fraction levels off 
or continues to increase. The efficiency will also probably 
depend on the structures within the gas disc. A more mas- 
sive disc, that is likely to excite lower order modes, will 
probably have a different effic iency to a lower mass disc 
that excites higher or der modes llLaughlin fc R6zvczkall996t 
lLodato fc Rice1l2004) . 

The above result shows that the particles that become 
strongly concentrated in the self-gravitating, gaseous, spiral 
structures, may achieve densities that could lead to plan- 
etesimal formation through direct gravitational collapse. To 
test that this is not simply a numerical effect, we carried out 
some additional simulations. We repeated the above simu- 
lation changing the gravitational softening for the solids to 
10 _1 au and 10~ 3 au. To a certain extent, softening acts 
to effectively smooth out the mass distribution. If the soft- 
ening is too large, this can act to suppress clumping. If it 
is too small, the 'graininess' of the mass distribution can 
exacerbate the gravitational dynamics, increasing the veloc- 
ity disperion of the particles. In our simulations, the par- 
ticle separation at the beginning of clump formation was 
~ 10~ 2 au, comparable to our original softening of 10 -2 
au. Although clumps formed for all of the softenings con- 
sidered, the maximum clump density was much lower when 
the largest softening was used, consistent with clump forma- 
tion being suppressed for large softenings. For the smallest 
softening considered, clump formation was delayed relative 
to the larger softenings. This delay in fragmentation is, as 
mentioned earlier, a consequence of the increased velocity 
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Figure 3. Fraction of 'planetesimal' particles located in clumps 
with densities exceeding 6 X 10 — 8 g cm -3 . The first sign of clump 
formation occurs about 20 years after the solid particles are intro- 
duced into the simulation with about 15% of the 'planetesimal' 
particles located in dense clumps after ~ 80 years. 

dispersion of the solid particles due to softening not be- 
ing large enough to effectively reduce the 'graininess' of the 
particle distribution. We also performed a simulation with 
1500 cm particles, also with an initial surface density 100 
times smaller than that of the disc gas, and another with 
150 cm particles, but with an initial surface density 1000 
times smaller than that of the disc gas. In both of these 
simulations the gravitational softening was 10 -2 au. 

Figure 2] shows the surface density structure of the sim- 
ulatio n with 1500 cm particles. As illustrated bv lRice et alJ 
(2004) particles of this size are largely decoupled from the 
disc gas, are essentially unaffected by gas drag, and hence do 
not become significantly concentrated in the gaseous spiral 
structures. The spiral structures in the 'planetesimal' disc 
effectively match that of the gas disc and are produced pri- 
marily by the gravitational potential of the gaseous disc, and 
not by gas drag. 

Figure |S] shows the surface density structure of the sim- 
ulation with 150 cm particles, with an initial surface density 
1000 times smaller than that of the gas disc. In this case, 
the particles do become concentrated in the gaseous spiral 
structures, producing very narrow spiral structures in the 
'planetesimal' disc. However, unlike the simulation with the 
higher initial surface density, the densities achieved do not, 
in this simulation with the lower initial surface density, lead 
to fragmentation and clump formation in the 'planetesimal' 
disc. It appears, therefore, that the gravitational collapse of 
solids requires particles that will become strongly concen- 
trated in the spiral density and that have a sufficiently high 
initial surface density, relative to the surface density in the 
gas disc. 



4 DISCUSSION AND CONCLUSIONS 

The simu l ations presented here extends the work of 
iRice et all i2004h which showed that particles of a partic- 
ular size (in their case 100 cm particles) , could become very 
strongly concentrated in the spiral structures present in a 
quasi-steady, self-gravitating gaseous protoplanetary disc. 

In this work we show that if the initial 'planetesimal' 
surface density is sufficiently high, the 'planetesimal' disc 
may undergo fragmentation to produce bound objects. It is 
this process that may provide a mechanism for producing 
kilometre-sized and larger planetesimals. In the simulations 
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Figure 4. Surface density structure of the 'planetesimal' disc 
consisting of 1500 cm particles ~ 100 years after the 'planetesi- 
mal' disc is introduced into the simulation. The scale shows the 
logarithm of the surface density, S p i, with the scale covering 
-0.6 < log(E/gcm" 2 ) < 2.4. 




Figure 5. Surface density structure of the 'planetesimal' disc 
consisting of 150 cm particles ~ 80 years after the 'planetesi- 
mal' disc is introduced into the simulation. The scale shows the 
logarithm of the surface density, E p i, with the scale covering 
-1.6 < log(S pl /gcm- 2 ) < 1.4. 

here each bound fragment consists of ~ 100 pseudo-particles 
and hence have masses of ~ 0.5 M cart h- These clumps are 
extremely strongly bound with the magnitude of the gravita- 
tional potential energy many times greater than the kinetic 
energy. It is therefore likely that each bound clump may rep- 
resent a number of smaller fragments, rather than a single 
bound object. 

We are therefore not suggesting that this process would 
lead, immediately, to earth mass-like objects, but simply 
that the densities and velocity dispersions of the particles 
concentrated in the self-gravitating spiral can lead to frag- 
mentation of the 'planetesimal' disc. 

We should stress, however, t hat the process here differs 
considerably from t he standard Goldreich fc Wardl (Il973l) 
mechanism (see also lYou din 2005) which requires a low ve- 
locity dispersion for fragmentation of the 'planetesimal' disc. 
In our simulations, fragmentation occurs when the velocity 
dispersion of the particles is comparable to the sound speed 
of the disc gas. It is the enhanced density that leads to the 
fragmentation, rather than the reduced velocity dispersion. 

It does appear, however, that for fragmentation to oc- 



cur, the surface density of particles in the appropriate size 
range must be relatively high. This could occur if, while the 
disc is still massive, particle growth up to ~ metre sizes 
occurs through collisional growth. Growth beyond metre 
sizes may, however, be diffi cult due to the stickin g efficiency 
becoming relatively small iSupulver et al.lll997f) and there 
would therefore be a pile up of these particles. Once the 
surface density of these particles reaches a critical value, 
subsequent growth would then occur through gravitational 
collapse triggered by the density of these particles being 
enhanced in the spiral structures of the gaseous disc. Ad- 
ditionally, the dependence of this process on the solid to 
gas ratio is consistent with planet for mation being more 
likely in higher metallic ity environments llSantos et al .120041 : 
iFischer fc Valentil2005l) . 

The process described here potentially solves a major 
problem in the standard planet formation scenario. Rather 
than rapidly migrating into the central star, ~ metre- 
sized particles become concentrated in self-gravitating spiral 
structures. The densities achieved can then lead to plan- 
etesimal formation via direct gravitational collapse. In this 
scenario, kilometre-sized planetesimals form very early, re- 
moving a major bottleneck in the planet formation process. 
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